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In the strong field molecular tunneling ionization theory of Tong et al. [Phys. Rev. A 66, 
033402 (2002)], the ionization rate depends on the asymptotic wavefunction of the molecular orbital 
from which the electron is removed. The orbital wavefunctions obtained from standard quantum 
chemistry packages in general are not good enough in the asymptotic region. Here we construct a 
one-electron model potential for several linear molecules using density functional theory (DFT). We 
show that the asymptotic wavefunction can be improved with an iteration method and after one 
iteration accurate asymptotic wavefunctions and structure parameters are determined. With the 
new parameters we examine the alignment-dependent tunneling ionization probabilities for several 
molecules and compare with other calculations and with recent measurements, including ionization 
from inner molecular orbitals. 

PACS numbers: 33.80.Rv, 42.50.Hz 
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I. INTRODUCTION 

Tunneling ionization of molecules in strong infrared 
fields is the first step in many interesting strong- 
field phenomena such as high-order harmonic generation 
(HHG), emission of high-energy above-threshold ioniza- 
tion (HATI) electrons and non-sequential double ioniza- 
tion (NSDI). Essential understanding to these processes 
is the angle- dependent ionization probability P{9) for a 
molecule fixed in space, where 9 is the angle between 
the molecular axis and the polarization direction of the 
laser's electric field. Since molecules are generally not 
fixed in space, i.e., not at a fixed alignment and/or ori- 
entation, experimental determination of P{9) from par- 
tially aligned molecules requires additional assumptions. 
Alnaser et al. first determined P{9) from NSDI pro- 
cesses where the alignment of the molecule is determined 
by Coulomb explosion of the molecular ions. P{9) can 
also be determined by ionizing partially aligned molecules 
d, 0] , or by measuring the angular distribution of elec- 
trons removed by a circularly polarized laser 0, In 
both methods the alignment of the molecular axis is de- 
termined by Coulomb explosion when the molecular ion 
is further ionized by an intense circularly polarized laser. 
In all of these measurements, the P{9) is not determined 
directly for a fixed angle and some approximations are 
used in order to determine the alignment-dependent ion- 
ization probability. 

Theoretically, P{9) can in principle be obtained di- 
rectly from numerical solution of the time-dependent 
Schrodinger equation (TDSE). However, even for the 
simplest H^, the P{9) obtained from solving TDSE 
by different groups still exhibits relatively large differ- 
ences. While calculations of P{9) for interesting multi- 
electron molecular systems have been carried out using 



the time-dependent density-functional theory (TDDFT) 
(see, for example, @), the accuracy of these calculations 
is difficult to evaluate. Furthermore, these calculations 
are rather time-consuming. Beside these ab initio ap- 
proaches, alignment-dependent tunneling ionization rate 
for molecules can be calculated using simple models such 
as the molecular strong field approximation (SEA) ^, ^Sj] , 
or the molecular tunneling ionization theory [9| . The lat- 
ter is the simplest and is a generalization of the tunnel- 
ing model of Ammosov, Delone and Krainov (ADK) [l3| 
for atoms. In the molecular tunneling ionization model 
(MO- ADK) of Tong et al. Q, the ionization rate for a 
molecule aligned at an angle 9 with respect to the laser 
polarization axis is given analytically. The ionization 
rate depends on the instantaneous electric field of the 
laser, the ionization potential of the molecule and some 
structure parameters of the orbital wavefunction in the 
asymptotic region. Subsequent further extension of the 
MO- ADK theory can be found in [uHIl- 

In Tong et al. @, the structure parameters are ex- 
tracted from molecular wavefunctions calculated using 
the multiple scattering method 14]. However, these 
days molecular wavefunctions are more easily accessible 
from quantum chemistry packages such as GAMESS [l^ , 
GAUSSIAN [ill and others. Thus it is desirable to ob- 
tain structure parameters from the asymptotic behavior 
of orbitals calculated from such packages. This was car- 
ried out for CO2 by Le et a l IItII and for other molecules 
by Kjeldsen and Madsen [18] . Unfortunately, molecu- 
lar orbitals from these chemistry packages are calculated 
using gaussian basis functions and they are not suitable 
for representing the exponential decay of the wavefunc- 
tion at large distances. As more accurate experimental 
data are becoming available, it is essential to redeter- 
mine these structure parameters more accurately. Since 
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the asymptotic wavefunction does not contribute much 
to the total energy of a molecule, one cannot efficiently 
improve the asymptotic wavefunctions by enlarging the 
size of the gaussian basis directly. 

In this paper, we describe how to improve the asymp- 
totic wavefunction where the structure parameters are 
extracted. Our input consists of wavefunctions of all 
the occupied orbitals obtained from GAMESS or GAUS- 
SIAN. We then construct a single-active-electron model 
potential and solve the time-mdependent Schrodinger 
equation to obtain the molecular orbital wavefunctions 
by an iterative procedure. The details of the method 
are given in Section II. We then apply the method to 
redetermine all the structure parameters previously pub- 
lished in 0, and adding structure parameters for some 
inner orbitals. We also determine the structure parame- 
ters for a number of systems that have been investigated 
experimentally. Using these new structure parameters we 
examined the alignment dependence of ionization proba- 
bilities for several systems. In most cases, the new results 
do not differ much from what were presented in Tong et 
al. [9]. However, there are differences in some molecules. 
The strong deviation in CO2 has been reported recently 

m. 



II. THEORETICAL METHODS 

The theory part is divided into three subsections. We 
first present the method of generating a single-active- 
electron model potential for linear molecules. We then 
discuss how to calculate the wavefunctions by solving the 
time-mdependent Schrodinger equation with B-spline ba- 
sis functions. We will also briefly describe how to extract 
the structure parameters in the MO-ADK theory. 



A. Construction of single-active-electron model 
potentials for linear molecules 



pressed in single-center expansion as 



V{r,e) = J2 Mr)Pi{cose). 



(1) 



(=0 



Here, w/(r) is the radial component of the model poten- 
tial and Pi{cos9) the Legendre polynomial. Typically we 
choose Imax — 40. The radial potential is given by 



v,ir) = vr'ir)+vfir)+vrir), 



(2) 



where the first two terms represent the electrostatic po- 
tential and the last term is the exchange interaction. 
The electron- nucleus interaction y^^'^ (r) can be written 

as 



(3) 



1=1 



where i runs over the Na atoms in the molecule. Without 
loss of generality, we assume that linear molecules are 
aligned along the z-axis, then W;(r) can be expressed as 



for Zi > 



vlir)=< ■> ^. ^ (4) 



with r^=min(r,|2;i|), r5^=max(r,|zi|). Here and Zi are 
the nuclear charge and the z coordinate of the ith atom, 
respectively. 



The partial Hartree potential vf (r) is given by 



21 + 1 



ai{r )r 



J+i 



(5) 



/o -> 
with r<=min(r,r'), r>=max(r,r'). Here a/(r') is 



2/ + 1 

aiir') = J ^p{r\e')Piicose')d{cos0'), (6) 



Single-active-electron model potential approach has 
been widely used for describing atoms in strong-field 
physics (see, for example, [20|). This approach has 
also been used for molecular targets recently [H, |2]| . 
The one-electron model potential consists of two parts: 
electrostatic and exchange-correlation terms. It is well- 
known that the traditional local-density approximation 
(LDA) for the exchange-correlation potential does not 
give the correct (— 1/r) potential in the asymptotic re- 
gion where the structure parameters are to be extracted. 
In this paper, we follow Abu-samha and Madsen [2l| and 
use the LB potential, proposed by Leeuwen and Baerends 
[221] . which will give the correct asymptotic —1/r behav- 
ior for neutral atoms and molecules. We note that a sim- 



ilar LB potential, called LBa [23[, has also been use d by 
Chu and collaborators in their TDDFT approach d, [2J . 
For linear molecules, the model potential can be ex- 



where p is the total electron density in the molecule and 



^V. (7) 



Here i runs over all the TVe electrons in the molecule. The 
wavefunction of each molecular orbital can be obtained 
from quantum chemistry packages such as GAMESS 
or GAUSSIAN [M]. 

For the partial exchange potential, it is written as 

27 + 1 /"i 

vr{r) = j ^ Ve.,air,9)Pi{cos9)dicos9), (8) 

where 
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Here V^J^J^ {r, 9) is the LDA potential for an electron with 
spin a 



1/3 



(10) 



where 



E 

i=l 



^ [ \^,,{r,e,^)fd^. (11) 



Here i runs over the N„ electrons that have the same 
spin as the active electron. The gradient correction term 
is given by [l^] 



l + 3l3xair, 0) smh-'ix^r, 0)) 



(12) 



where Xcr(?', ^) — \'^Pa-{r,9)\pa'^^^{r,0). The parame- 
ters a and /? are chosen to be 1.0 and 0.05, respectively 
throughout this paper. We note that for more accurate 
binding energies, the correlation potential should be in- 
cluded into Eq. (9). In the so-called LBa model, the two 
parameters a and /3 are usually chosen to be 1.19 and 
0.01, respectively (see, [23|). 



B. Calculation of molecular wavefunctions by 
solving the time-independent Schrodinger equation 

With the model potential constructed in the previous 
subsection, the wavefunction for the active electron in a 
linear molecule can be obtained by solving the following 
time-mdependent Schrodinger equation 

77,zVi"^(r) EE [-iv2 + V{r,0M^rHr) = i?„^i™)(r) 

(13) 

where ipn"'^ and £'4™'* are the eigenfunction and eigen- 
value, respectively. 

Using single-center expansion for the electronic wave- 
function 



e^(r)=E^ll™(^,^) 



(14) 



1=0 



where Yim{0T^) are the spherical harmonics, the radial 
wavefunction can be constructed with B-splines [1^ 



(15) 



4=1 



Substituting Eqs. (1), (14) and (15) into Eq. (13) and 
then projecting onto the BiYj*^ basis, we obtain the fol- 
lowing matrix equation 



HC = ESC 



(16) 



where 



il.i'l' 



f B,{r)Yi*„,{0,v)Hei 
"'0 Ja 

Bi'Yiim{0, (p)dr snY0d0dLp 



(17) 



SiiA'v = 



B,{r)B,,{r)dr (18) 



The eigenfunctions and eigenvalues are obtained by di- 
agonalizing Eq. (16). 



C. Extracting asymptotic structure parameters 

In the asymptotic region, typically only a few terms 
in the single-center expansion Eq. (14) are important. 
Following Tong et al. @ , we write the wavefunction of a 
linear molecule as 



(19) 



In the MO-ADK theory [9|, the radial functions in the 
asymptotic region are fitted to the following form 



where is the asymptotic charge, k, = 
ionization energy. 



(20) 



/2/p, Ip is the 



III. RESULTS AND DISCUSSION 

A. On the quality of the model potential and the 
iteration procedure 

In this paper, the single-active-electron model poten- 
tial [see Sec. II A] is created with the DFT, in which 
the exchange potential is constructed with the exchange- 
only LDA potential and the LB model potential (or 
LDA-I-LB). First we check the quality of this model po- 
tential if the molecular orbitals obtained from the stan- 
dard quantum chemistry package GAMESS pj| are used 
as the input. 

In Fig. 1(a), we compare the present r- weighted LB 
potential with the empirical model potential of Tong et 
al. [131 for Ne. For clarity we plot the effective charge, 
defined as r x V{r). The two potentials agree well in 
the small r region. However, there are significant differ- 
ences at large r. For neutral atoms, the effective charge 
should approach -1 at large r. If the LB potential is cal- 
culated directly using the molecular wavefunctions from 
GAMESS (dashed line) the effective charge exhibits os- 
cillations and then drops rapidly with r. This undesir- 
able behavior is due to the incorrect electron density, 
which in turn is due to the limitation of the gaussian 
basis, calculated from GAMESS in the large r region. 
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FIG. 1: (Color online) (a) Effective charge of Ne with and 
without the iteration (see text). Model potential is from [20l |. 
(b) Effective charge of CO2 along the molecular axis. 



To correct this error, we perform one more iteration on 
the potential: Firstly, an initial model potential is gen- 
erated using the Hartree-Fock (HF) wavefunctions ob- 
tained from GAMESS. From this initial potential, more 
accurate wavefunctions are obtained by solving Eq. (13) 
with B-spline basis. Then, a new model potential is con- 
structed from these new wavefunctions. From Fig. 1(a), 
we observe that the effective charge obtained after one 
iteration (solid line) shows the correct asymptotic behav- 
ior. The same procedure can be applied to molecules. In 
Fig. 1(b), we show the model potential of CO2 along the 
molecular axis, with and without one iteration. It con- 
firms that the asymptotic behavior of the model potential 
is correct after one iteration. We comment that in the 
case of CO2 diffuse functions have been included in the 
basis sets. Clearly, this alone is insufficient for obtaining 
accurate electron density (or potential) at large r. 



B. Extracting molecular structure parameters for 
the MO-ADK theory 

Once the model potential is obtained, the eigen- 
function and eigenvalue can be calculated from solv- 
ing Eq. (13). In Table I, binding energies of rare gas 
atoms obtained using the present method are compared 
to those from Ref. 26] and the experimental values. Our 



TABLE I: Comparison of calculated ionization energies of rare 
gas atoms in the exchange-only LDA+LB model and experi- 
mental values. 



Atom 


LDA-hLB ( 


a.u.) 


Ip (a.u.) 


He 


0.786 




0.904" 




0.796* 






Ne 


0.722 




0.793" 




0.725'' 






Ar 


0.524 




0.579" 




0.528'' 






Kr 


0.499 




0.515" 


Xe 


0.469 




0.446" 


"Reference 


Si 






''Reference 


3 






TABLE II: 


Equilibrium distances, ionization energies calcu- 


lated in the exchange-only LDA-I-LB model, and 


experimental 


vertical ionization potentials for several linear molecules. 


Molecule 


R(A) 


LDA+LB (eV) 


Ip (eV) 




1.058 


29.99 


29.99 


D2 


0.742 


13.65 


15.47 


N2 


1.098 


14.99 


15.58 


O2 


1.208 


10.62 


12.03 


F2 


1.412 


16.03 


15.70 


S2 


1.889 


10.36 


9.36 


CO 


1.128 


13.22 


14.01 


NO 


1.151 


9.14 


9.26 


SO 


1.481 


9.37 


10.29 


CO2 


1.163 


14.63 


13.78 


C2H2 


1.203 (Rcc) 


11.19 


11.41 




1.058 (Rch) 






HF 


0.917 


15.03 


15.77 


HCl 


1.275 


11.41 


12.75 


HCN 


1.067 (Rch) 


13.46 


13.80 




1.159 (Ron) 







method uses the same approximate exchange potential 
as in Ref. [26]. The two calculations agree in general, 
but small discrepancies do exist with experimental val- 
ues. These discrepancies can be reduced if correlation 
potential is included in Eq. (9). This fact has been well 
documented in Ref. [2^. 

In Table II, we compare the ionization energies from 
the present calculations with experimental vertical ion- 
ization energies for several linear molecules. The equi- 
librium distances of these molecules are also listed. The 
agreement between the calculated and experimental val- 
ues are good. Again we comment that the exchange- 
only LDA-t-LB potential are used in our calculations. For 
hi gher precision, correlation potential should be included 

Slllillii. 
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TABLE III: The newly fitted Cim coefficients vs values from 
earlier references, P. [itI. [29l. [30| . 



Molecule 




^\ra 














■p"+ ^ 


/I KO 

4.0Z 




A f?0 




A AO 

U.Uo 










4.0 ^ 




A Att 




A AA 
U.UU 






roi 

m 




1. (^O 




All 
U.ii 




A AA 
U.UU 










OKI 

z.oi 




A AiS 




A AA 
U.UU 






FaI 

m 




11c: 

i.io 




A A^l'T' 

U.Ud / 




A AA1 
U.UUl 






ronl 




Z.DO 




1 1 A 
i.iU 




A Afi 

U.Ud 










no 




A '7'C 

U. / 




A A/1 

U.U4 






roi 

m 


02(71-9) 






A CO 

U.oz 




A AO 

U.Uo 














n £io 




A AO 

U.Uo 






FaI 


"PT f''n- A 






l.Zl 




U.lo 














1.17 




A 1 

U.lo 








b2(7ra) 






l.OY 




A 1 '7 
U.l / 














U.ol 




A PiT 
U.U 1 






roi 




QO 

z.oz 


1 i^O 
i.Dz 


A 00 
U.OZ 


U.i 1 


A AC 

U.Uo 










1 /I Q 


U. / D 


U.zo 


n no 


U.UU 












A 01 

U.zi 


A 00 

U.OO 


U.Oz 


A AO 

U.Uz 












A 00 

U.zz 


A -1 1 

U.4i 


U.Ui 


A AA 
U.UU 






rnl 

m 


oU(7r ) 




A 00 


A T1 

yj.i L 


0.05 


A AC 

U.Uo 












A /I 1 


A Ql 

-U.ol 


U.Ui 


A AA 
U.UU 






rnl 
19] 


C02(7rg) 






1.97 




0.40 




0.04 










2.88 




1.71 




0.43 


m 


C2H2(7r„) 




1.16 




0.18 




0.02 










1.14 




0.27 




0.04 




[30] 


HF(7r) 




0.88 


0.03 


0.02 


0.01 








HCl(7r) 




1.23 


0.01 


0.05 


0.01 


0.01 






HCN(7r) 




1.50 


0.09 


0.24 


0.02 


0.02 







With the new wavefunctions, we re-evaluate the struc- 
ture parameters for a number of linear molecules. Table 
III lists the newly fitted C/„i coefficients with those listed 
in Tong et al. |9| and in others, if available. These pa- 
rameters will be used to obtain the alignment-dependent 
tunneling ionization rates, following the MO-ADK theory 

i- 

C. Comparison of alignment-dependent ionization 
probabilities between MO-ADK and other 
calculations 

Using the improved structure parameters tabulated in 
Table III, we now use the analytical formula in Tong et al. 
to obtain alignment-dependent tunneling ionization 
probabilities for selected molecules that have also been 
carried out by other methods. The results are shown in 
Fig. 2. For simplicity, all the probabilities are normalized 
to 1.0 at the peak. First, we comment that for N2, O2, 
F2 the normalized probabilities obtained using the new 
structure parameters do not show noticeable differences 
compared to the probabilities calculated using old struc- 



ture parameters. From Table III, we note that the struc- 
ture parameters for these three molecules do not change 
much. We emphasize that in calculating MO-ADK rates, 
one should always use the experimental vertical ioniza- 
tion energy since the tunneling ionization rate depends 
exponentially on the ionization potential. In Figs. 2(d) 
and 2(e), we notice that, interestingly, the MO-ADK re- 
sults using the new Cim give stronger angular dependence 
than the old ones for both H J and II2 . This is the result 
of the relatively larger C2m as compared to Com in the 
present calculations. For H^, the present result lies be- 
tween the two calculations from solving TDSE. For H2, 
we compare the new results with those from SFA, and 
the two agree quite well. For C2H2, the new MO-ADK 
result agrees with the SFA, but differs from the older 
MO-ADK [30]. We comment that in the SFA calcula- 
tion, wavefunctions directly from the GAMESS code are 
used. In general, SFA calculations yield incorrect total 
ionization rates. Empirically, however, the normalized 
alignment dependence from the SFA appears to be in 
agreement with the present MO-ADK. In presenting the 
SFA results, we always use the renormalized ones. We 
further comment that in SFA and other ab initio calcu- 
lations, ionization probability or rate for each alignment 
angle is calculated independently. In the MO-ADK the- 
ory, the alignment dependence is obtained analytically 
after the structure parameters are obtained. 

In recent years, ab initio calculations of molecular ion- 
ization by intense lasers have been carried out by solving 
the TDDFT [6,,^, |31|]. These calculations include all 
the electrons in the molecule. Comparing to MO-ADK, 
in general, these calculations tend to give larger prob- 
abilities at angles where the ionization is small, see N2 
near 90° and O2 and F2 at angles near 0° and 90°. For 
C2II2, on the other hand, the TDDFT result is smaller 
at smaller angles than the present one. For this system, 
it was carried out by a different group 'si']. Based on 
these results we can say that the alignment dependence 
of the ionization probabilities obtained from MO-ADK 
and from TDDFT are in reasonable agreement. However, 
we mention that probabilities in Fig. 2 from MO-ADK 
include ionization from the HOMO only, while the many- 
electron TDDFT calculations show significant contribu- 
tions from the inner orbitals. More on the comparison 
between MO-ADK and TDDFT will be given later. 



D. Alignment dependence of ionization rates from 
HOMO, HOMO-1 and HOMO-2 orbitals 

Recently, strong field ionization phenomena involving 
inner orbitals of molecules have been reported widely 
[s^ - iSTj . This is somewhat surprising since tunnel- 
ing ionization rate decreases very rapidly with the in- 
crease of ionization potential. However, molecular tun- 
neling ionization rates depend on the symmetry of the 
orbital wavefunctions. For alignment angles where P{0) 
is near the minimum for the HOMO but where HOMO-1 
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FIG. 2: (Color online) Normalized alignment dependence of ionization probability, (a) N2 at laser intensity of 10^* W/cm'^; (b) 
O2 at lO^"* W/cm^; (c) F2 at 2 x 10^'' W/cm^ (d) at 5 x lO" W/cra-Je) H2 at 2.3 x lO^"* W/cm^; (f) C2H2 at 5 x lO" 
W/cm^ TDDFT" from Telnov et al. TDDFT*' from Otobe et al. TDSE" from Kamta et al. [H, TDSE*' from 

Kjeldsen et al. :33], Tong et al. from 0] and Lin et al. from [30l ]. 



is near the maximum, there is a good possibility that ion- 
ization from HOMO-1 can become comparable or higher 
than from HOMO. Indeed, contribution from HOMO- 
1 to high-order harmonic generation (HHG) from Nj 
molecules has been reported by McFarland et al. [SJ] 
when the molecules are aligned perpendicular to the po- 
larization of the probe laser. Le et al. [35| have suc- 
cessfully reproduced the experimental results by includ- 
ing HHG from HOMO and HOMO-1. Since tunneling 
ionization is the first step for all rescattering processes 
[38l - l4l| , including HHG 'w] , it is pertinent to investigate 



V{6) from inner orbitals as well. 

In Table IV, the binding energies of HOMO, HOMO-1 
and HOMO-2 for several molecules are shown. These en- 
ergies are compared to calculations using the LBa model 
and experimental values, to check the relative accuracy of 
the model we have used. We emphasize again that accu- 
rate experimental ionization energies, not the theoretical 
values in the Table, are used in calculating the MO-ADK 
rates. The extracted Cim parameters are given in Table 
V. Using these parameters and experimental ionization 
energies, the alignment dependence of ionization rates 
from different orbitals at a given peak laser intensity can 
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TABLE IV: Comparison of calculated binding energies of 
HOMO, HOMO-1 and HOMO-2 of Na, O2 and CO2 in the 
present exchange-only LDA+LB model. Those from the LBa 
model and experimental vertical ionization potential are also 
given. Energies are in electron volts. For CO, HCl and C2H2, 
only the energies of HOMO and HOMO-1 are considered. 



Molecule 




Spin orbital 


LDA+LB 


LBq 


Ip 


N2 




3cr9(HOMO) 


15.0 


15.5" 


15.6'' 




1 


7r„(H0M0-l) 


16.5 


16.9" 


17.2'' 




2cr„(HOMO-2) 


17.8 


18.5" 


18.7'' 


O2 




l7r9(HOMO) 


10.6 


12.8" 


12.3"= 




1 


7r„(H0M0-l) 


17.3 


17.4" 


16.7^= 




So-g (HOMO-2) 


17.1 


18.3" 


18. 2^= 


CO2 




l7r9(HOMO) 


14.6 


13.9^* 


13.8" 




l7r„(H0M0-l) 


18.3 


17.5"^ 


17.6'= 




3cr„(HOMO-2) 


16.8 


17.2'* 


18.1'= 


CO 




5cr(HOMO) 


13.2 




14.0'= 






l7r(HOMO-l) 


16.6 




16.9'= 


HCl 




27r(HOMO) 


11.4 




12.8^ 




5a(HOMO-l) 


15.0 




16.3-'' 


C2H2 




l7r„(H0M0) 


11.2 




11.4'= 




So-g (HOMO-1) 


15.7 




16.4'= 


"Reference 
''Reference 
'^Reference 
Reference 
'^Reference 
■'^Reference 


2^ 
45 
4S 

44 
4£ 





be readily calculated. 

In Fig. 3, we compare the ionization rates from N2, 
O2 and CO2 molecules, for the HOMO, HOMO-1 and 
HOMO-2 orbitals, at peak intensities indicated in the 
figure. Note that the angular dependence, P{d), reflects 
the symmetry of the molecular orbital quite accurately. 
Thus a a orbital tends to have the peak at 0° and a 
minimum at 90°, a Wg orbital has the peak near 45° and 
minimum at 0° and 90° , and a 7r„ orbital has a peak near 
90° and minimum near 0° (Deviations do occur, see the 
HOMO-1 of CO2 in Fig. 3(c)). These general behav- 
iors of ionization rates explain why HOMO-2 is bigger 
than HOMO-1 at small angles for N2, O2, and CO2, and 
why HOMO-1 is more important than HOMO at small 
angles for C2H2. Note that the relative ionization rates 
depend on laser intensities. The relative ionization rates 
for inner orbitals increases faster with increasing laser in- 
tensities. Using the parameters in Table V, their relative 
rates can be easily calculated using the MO-ADK model. 
We have also calculated the ionization rates using the 
molecular SFA. The relative alignment dependence from 
SFA in general agrees with those shown in Fig. 3. This 
is consistent with the findings in Le et al. [ssj . 

Fig. 4 shows the HOMO and HOMO-1 ionization rates 
for asymmetric diatomic molecules CO and HCl. There 
are recent experiments and other theoretical calculations 



TABLE V: The Ci coefficients of HOMO, HOMO-1 and 
HOMO-2 for N2, O2 and CO2 and of HOMO, HOMO-1 for 
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available for these two molecules [sB] • For both sys- 
tems, the predictions from MO-ADK are also compared 
to results from SFA. Refer to Table IV, we note that 
the difference in binding energies between HOMO and 
HOMO-1 in CO is 2.9 eV, and 3.5 eV for HCl. First 
we examine the 6'-dependence predicted by MO-ADK in 
Figs. 4(a) and 4(b). The HOMO of CO is a cr orbital, its 
P{9) drops rapidly from 0° to 90° and stays relative flat 
at larger angles. The HOMO-1 is a tt orbital and its P{9) 
peaks near 90°. For HCl, the HOMO is a tt orbital and it 
peaks near 90°. For the HOMO-1, it is a a orbital and its 
P{9) drops steadily till near 90°. Interestingly, its P{9) 
increases rapidly from 90° to 180°, making it almost like 
a symmetric molecule. 

Why are the a orbitals of the two molecules so differ- 
ent? It is due to the degree of asymmetry in the wave- 
functions. Such asymmetry is reflected in the C; coeffl- 
cients in Table V. For CO (HCl), the first three coeffi- 
cients are 2.32, 1.62, 0.82 (0.10, 2.64, 0.57) for ^=0, 1 and 
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FIG. 3: (Color online) Alignment dependence of ionization rates of HOMO, HOMO-1 and HOMO-2 for N2, O2 and CO2 and 
of HOMO and HOMO-1 for C2H2. (a) N2 at laser intensity of 1.5 x 10^* W/cm^; (b) O2 at 1.3 x lO" W/cm^; (c) CO2 at 
1.1 X 10^" W/cm^; (d) C2H2 at 1.5 x lO^"* W/cm^ 



2, respectively. For HCl, there is one dominant 1=1 com- 
ponent only, thus the ionization rate is nearly symmetric. 
For CO, the two coefficients for 1=0 and 1 are compara- 
ble, the wavef unction along the axis for 6 = and 6 = it 
has the ratio (2.32-|-1.62)/(2.32-1.62)=5.6. This gives a 
ionization rate ratio of 32, close to the value 50 read off 
from Fig. 4(a). 

In Figs. 4(a) and 4(b), the 6'-dependence from SFA is 
different from MO-ADK. Reeah that in MO-ADK, static 
ionization rate was calculated, thus a molecule is AB or 
BA with respect to the fixed electric field will have differ- 
ent rates. For a linearly polarized laser pulse, the direc- 
tion of the electric field changes after each half cycle, thus 
the cycle-averaged rates for AB and BA are identical. To 
compare the SFA rate with the MO-ADK rate at an angle 
6, we have to average the rates from the latter at 6 and 
ir-d. These "symmetrized" ionization rates are denoted 
by MO-ADK-S in Fig. 4. By comparing the rates from 
SFA and MO-ADK-S, we found in Fig. 4(c) that the two 
models agree well for CO. For HCl, the relative rates for 
HOMO-1, normalized to HOMO, are about a factor of 
two larger from SFA than from MO-ADK. We comment 
that if ionization is measured using circularly polarized 
light, the static MO-ADK rate can be compared directly 



with the rate calculated using SFA. 

The results of Figs. 3 and 4 show that at alignment an- 
gles where tunneling ionization from the HOMO is large, 
contributions from HOMO-1 or other inner orbitals are 
negligible. At alignment angles where HOMO is near the 
minimum, if the HOMO-1 (or even HOMO-2) is near the 
maximum, then these inner orbitals may become impor- 
tant. Since the relative tunneling ionization rates also 
depend on the peak laser intensity, when multiple or- 
bitals contribute to strong field phenomena, the inten- 
sity dependence may become prominent. Experimen- 
tally, such multiple orbital effects have been observed in 
HHG from N2 when molecules are alig ned perpendicu- 
lar to the laser's polarization axis |34l. Issj . The inner 
orbitals have been shown to become important in HHG 
from CO2 when the molecules are aligned parallel to the 
laser's polarization axis [s^l- Comparing to single photon 
ionization, strong field ionization tends to be more selec- 
tive by ionizing the HOMO. For single photon ionization, 
cross sections for HOMO, HOMO-1 and HOMO-2 in gen- 
eral have comparable values and often cross sections from 
inner orbitals are higher, see e.g., [i^ for CO2. 

There are few theoretical alignment-dependent ioniza- 
tion rates from inner orbitals available to compare with 
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FIG. 4: (Color online) Alignment dependence of ionization rates of HOMO and HOMO-1. (a) CO at laser intensity of 4 x 10^'' 
W/cm^; (b) HCl at 2 x 10^'' W/cm^ (c) CO at 4 x lO^"* W/cm^; (d) HCl at 2 x lO" W/cm^ MO-ADK-S is the averaged 
MO-ADK rate for angles 6 and n -9. 



the predictions of the MO-ADK theory presented here. 
There is an exception, however, CO2. In [47 ] ionization 
rates from HOMO, HOMO-1 and HOMO-2 have been 
calculated starting from the multielectron perspective. 
In Fig. 5(a) we compare the rates from MO-ADK with 
those from j47| at the uncoupled channel approximation. 
The two sets of calculations are normalized at the peak 
of the HOMO curve. We note that the 6'-dependence 
agrees well for each orbital. For the HOMO, the agree- 
ment is "perfect" . The rates for the inner orbitals are 
larger from [l^ than from MO-ADK. Part of the reason 
of the larger difference in the HOMO-1 rate could be due 
to the difference in the ionization energy used. In (47| . 
the energy difference between HOMO-1 and HOMO was 
taken to be 3.53 eV, while in MO-ADK, the difference 
was taken to be 3.80 eV from the experimental values in 
Table IV. For the HOMO-2 the energy used is the same 
for the two calculations. The alignment dependence of 
ionization rates for the three orbitals have also been cal- 
culated in [131 and the comparison with the present MO- 
ADK is given in Fig. 5(b), again by normalizing at the 
peak value of the HOMO. In this case the differences are 



larger. In [s?!, the ionization rates were calculated using 
Coulomb corrected SFA plus sub-cycle dynamics. The 
TDDFT method has also been used to obtain ionization 
probabilities from different orbitals [1, [2^ . For CO2 , the 
predicted alignment dependence for the HOMO, HOMO- 
1 and HOMO-2 as shown in Fig. 3 of @ do agree with the 
present Fig. 3(c), including that the peak for HOMO-1 is 
not at 90° . However, we should comment that in N2 and 
O2, the ali gnm ent dependence using the same TDDFT 
method in f2»| does not agree with Figs. 3(a) and 3(b) 
shown for these two molecules. 



E. Comparisons with Experiments 

Fig. 6 shows the normalized alignment dependence of 
ionization probability of N2, O2, H2 and HCl. From 
Figs. 6(a) and 6(b), we comment that the normalized ion- 
ization probability of N2 calculated from MO-ADK the- 
ory using the old and the newly fitted coefficients agree 
quite well (no visible difference in the plot). Compared to 
the experiment of Pavicic et al. , the MO-ADK theory 
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FIG. 5: (Color online) Comparison of ionization rates of HOMO, HOMO-1 and HOMO-2 of CO2 at peak laser intensity of 
1.5 X 10^* W/cm^. The solid lines are from MO-ADK and the dashed lines are from Spanner and Patchkovskii [i^] (a) and 
from Smirnova et al. 37] (b). 



shows some differences. But the difference is considered 
acceptable. Note that the determination of alignment de- 
pendence from the experiment has angular average which 
was not included in the theory curve. Take the experi- 
mental result as reference, the TDDFT result (see Fig. 2) 
is better than the MO-ADK for N2. For O2, it is the 
other way around. The same comparison for CO2 has 
been addressed in an earlier paper 19]. In that case, the 
old MO-ADK results were found to be inaccurate due to 
the inaccuracy of the old C; parameters. In [19] it was 
further concluded that the experimental P(0) from [3(] ap- 
pears to be too narrowly peaked. We note that the new 
result from [l^ also does not agree with the experiment. 
However the authors suspect that the discrepancy is due 
to intermediate excitation channels were not included in 
their calculation. We tend to think that additional ex- 
periments are needed to help resolving this discrepancy. 

In Fig. 2(e) we show that the MO-ADK probabilities 
for H2 using the new structure parameters are differ- 
ent from using the earlier ones [9|. The new MO-ADK 
probabilities and molecular SFA agree well, see Fig. 6(c). 
Comparing to experimental data of Staudte et al. , the 
agreement is good in view that the theory curve has not 
included average over angular resolution. In [3| , the ratio 
of ionization rate for molecules aligned parallel vs per- 
pendicular, with respect to the polarization axis, were 
also determined at four intensities from 2 to 4.5x10^^ 
W/cm^ (for circularly polarized laser). The ratio from 
the present SFA (not shown) agree with the SFA model 
in that paper, and with the new MO-ADK ratio of 1.45 
(the old MO-ADK gives 1.15). We expect the theoretical 



ratio be reduced somewhat if angular average is incor- 
porated. We mention that a similar measurement at one 
intensity for laser wavelength of 1850 nm was reported in 
psj . which gives a ratio of 1.15. Interestingly, this ratio 
was reported to be 3.0 [i^ in another recent experiment, 
while the theory presented in the same paper gives a ratio 
of 2.1. We comment that the ratio is taken at the maxi- 
mum with respect to the minimum and thus sensitive to 
the angular average. Comparison of the rates over the 
whole angular range would be preferable. 



In Fig. 6(d), the P{e) of the HOMO-1 orbital in HCl 
reported in Ref. Q using circularly polarized light at the 
intensity of 1.4x10^* W/cm^ is shown. We compare the 
HOMO-1 result from the MO-ADK theory using the laser 
parameters given in the experiment, and by normalizing 
the data at 6'=0°. In Ref. Q, the alignment dependence 
for HOMO and HOMO-1 has also been reported using 
the TDDFT. The alignment dependence between MO- 
ADK and TDDFT calculations are quite similar, but our 
relative HOMO-1 probability is about a factor of three 
higher at the same laser intensity. The ionization proba- 
bility from both calculations drop much faster from 0° to 
90° when compared to the experiment. By introducing 
a small fraction of the contribution from the HOMO in 
the manner suggested in [5|, the MO-ADK theory can 
achieve a reasonable agreement with the experimental 
data from 0° to 90°, see Fig. 6(d). On the other hand, 
the agreement at angles larger than 90° is still not as 
good. 
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FIG. 6: (Color online) Normalized alignment dependence of ionization probability, (a) N2 at laser intensity of 1.5 x lO'^* 
W/cm^; (b) O2 at 1.3 x lO" W/cm^; (c) H2 at 2.3 x 10^* W/cm^; (d) HCl at 1.4 x lO" W/cm^ Linearly polarized lights for 
(a) and (b); Circularly polarized lights for (c) and (d). Exp.° from Pavicic et al. [3]; Exp.*" from Staudte et al. U and Exp.'^ 
from Akagi et al. 3- Additional symbols for (d), see text. 



F. Ionization probability of 

The ionization probability of has been calculated 
from solving the TDSE by different groups [H, HE [131 • 
It is of interest to compare the predictions based on MO- 
ADK with those from solving the TDSE. In Fig. 7, the 
normalized alignment-dependent ionization probability 
from the first four molecular orbitals of at the equi- 
librium distance are shown. The data for IscTg have been 
discussed earlier [loj . For ionization from lscr„, the two 
TDSE calculations and the MO-ADK agree quite well. 
For 2p7r„, the MO-ADK theory tends to peak at 90° while 
the TDSE result gives a peak closer to about 60°. For 
2p7rg state, the MO-ADK predicts a peak near 45° while 
TDSE calculation gives a peak at about 55°. Note differ- 
ent peak laser intensities are used for the ionization from 
each orbital. 

In Fig. 8(a) we show the dependence of normalized 
ionization probabilities vs the internuclear separation for 
the IsCTg.M states of with the molecular axis parallel 
to the polarization axis. The results are compared to the 
TDSE calculations of [5l|. By normalizing the proba- 
bility at i? = 2 a.u. for the 1s(t„, we find that there is 



a general good agreement between the TDSE result and 
from the MO-ADK. For the Isct^, the two calculations 
are normalized at R = 4.0 a.u.. For both calculations, 
the probabilities at R less than 3.5 a.u. are significantly 
smaller than at i? = 4.0 a.u.. In Figs. 8(b) and 8(c), 
the normalized alignment-dependent ionization rates are 
shown for different R. Clearly as R increases, the angular 
dependence becomes sharper. This is easily understood 
for a orbitals since the molecular orbital becomes more 
elongated along the molecular axis as R increases. The 
Ci coefficients are tabulated in Table VI to reflect how 
these parameters vary as R increases. 



IV. CONCLUSIONS 

In this paper we proposed a new method to obtain 
accurate molecular wavefunctions in the asymptotic 
region starting with molecular orbitals obtained from 
the widely used quantum chemistry packages such as 
GAMESS and GAUSSIAN. From these wavefunctions, 
the structure parameters in the molecular tunneling 
ionization theory (MO-ADK) of Tong et al. @ can be 
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FIG. 7: (Color online) Normalized alignment dependence of ionization probability of Hj. (a) Isag at laser intensity of 5 x 10^'* 
W/cm^; (b) lsa„ at 10^* W/cm^; (c) 2p7r„ at 10^^ W/cm^; (d) 2p7Vg at lO" W/cm^ TDSE" from Kamta et al. TDSE*' 
from Kjeldsen et al. ^ and TDSE'= from Telnov et al. [13]. 



accurately determined. Using these structure parame- 
ters, we re-examined the alignment-dependent tunneling 
ionization probabilities for a number of molecules, in- 
cluding ionization from HOMO-1 and HOMO-2 orbitals. 
The calculated tunneling ionization probabilities are 
compared to probabilities determined from experiments, 
and to several other more elaborated calculations. 
Since tunneling ionization is the first step for strong 
field phenomena involving molecular targets, these 
structure parameters are useful and thus are tabulated. 
The procedure for obtaining the structure parameters 
discussed in this paper is generally applicable to any 
linear molecules. Despite of its fundamental importance, 
accurate strong field alignment-dependent ionization 
probabilities are still not widely available. Experimental 
measurements as well as more advanced calculations 
tend to deal with different molecules and under different 



conditions, thus it is difficult to benchmark the accuracy 
of the theoretical models. While MO-ADK model is 
the simplest model for obtaining tunneling ionization 
rates, it appears that its predictions so far are in good 
agreement with most of the experimental data and with 
most the elaborate theoretical calculations. 
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